

explim=10;
for precond=0:1
    for n=2.^(1:explim)
		precision='double';
		if n > 2^(explim-1)
			precision='single';
		end
        L=3*n-2; %n+2n-1-1
        filename=sprintf('filters%sH2Dfilter%dprecond%dprecision%s.dat',filesep,n,precond,precision);
        if exist(filename,'file')
            H=load2DHFilter(filename,precision);
        else
            H=save2DHFilter(filename,n,precond,precision);
        end
        if 0
            Hhat=fftn(H.filter,[L L]);
            e0=zeros(n);
            T=zeros(n^2);
            for i = 1:n^2
                    e1=e0;
                    e1(i)=1;
                    e1=applyHFilter2D(e1,Hhat);
                    T(i,:)=e1(:); %symmetric operator
            end
            [u s v1] = svd(T);
            sd=diag(s);
            sbig=sd(1);
            ssmall=sd(n^2);
            svdConditionNumber=sbig/ssmall;
            [v2 d] = eig(T);
            dd=diag(d);
            dbig=dd(1);
            dsmall=dd(n^2);
            eigConditionNumber=dbig/dsmall;

            v=zeros(n^2,1);
            rr=0;
            rrm1=0;
            for i=1:1e4
                a=pi*rand(n^2-1,1);
                a(n^2-1)=2*a(n^2-1);
                for j=1:n^2
                    if j==n^2
                        r=1;
                    else
                        r=cos(a(j));
                    end
                    for k=1:j-1
                        r=r*sin(a(k));
                    end
                    v(j)=r;
                end
                rr=max(rr,norm(T*v));
                rrm1=max(rrm1,norm(T\v));
            end
            defConditionNumber=rr*rrm1;
        end
        tic
        powerConditionNumber = calculateConditionNumber(H,precision);
        toc
        %fprintf('n = %d, svdConditionNumber = %4.6f, powerConditionNumber = %4.6f, eigConditionNumber = %4.6f defConditionNumber = %4.6f\n'...
        %    ,n,svdConditionNumber,powerConditionNumber,eigConditionNumber,defConditionNumber);
        fprintf('n = %d, precond = %d,  powerConditionNumber = %4.6f\n'...
            ,n,precond,powerConditionNumber);
        %reportTestResult(n,abs((svdConditionNumber-powerConditionNumber))/svdConditionNumber,t1/t2,precision);    
    end
end